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Abstract 

A quantitative understanding of organism-level behavior requires pre¬ 
dictive models that can capture the richness of behavioral phenotypes, yet 
are simple enough to connect with underlying mechanistic processes. Here 
we investigate the motile behavior of nematodes at the level of their trans¬ 
lational motion on surfaces driven by undulatory propulsion. We broadly 
sample the nematode behavioral repertoire by measuring motile trajecto¬ 
ries of the canonical lab strain C. elegans N2 as well as wild strains and 
distant species. We focus on trajectory dynamics over timescales spanning 
the transition from ballistic to diffusive movement and find that salient 
features of the motility statistics are captured by a random walk model 
with independent dynamics in the speed and bearing. We show that the 
model parameters vary among species in a correlated, low-dimensional 
manner suggestive of a common mode of behavioral control and a trade¬ 
off between exploration and exploitation. The distribution of phenotypes 
along this primary mode of variation reveals that not only the mean but 
also the variance varies considerably across strains, suggesting that these 
clonal populations employ contrasting ’bet-hedging’ strategies for forag¬ 
ing. 
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Abbreviations: MSD, mean-squared displacement; VACF, velocity autocor¬ 
relation function; ACF, autocorrelation function; SDS, sodium dodecyl sulfate 


Introduction 

A ubiquitous feature of biological motility is the combination of highly stereo¬ 
typed movements in seemingly random sequences. Capturing the essential char¬ 
acteristics of motion thus requires a statistical description, in close analogy to 
the random-walk formulation of Brownian motion in physics. Indeed, the the¬ 
ory of random walks has been applied broadly to study biological motility, and 
has uncovered fundamental differences between physical Brownian motion and 
biological random walkers. For example, over long times, the random motility 
of many biological organisms appear to follow non-Brownian diffusive dynam¬ 
ics known as Levy flights - random walks in which the length distribution of 
uninterrupted flight segments follow a power law (as opposed to the exponen¬ 
tial distribution of Brownian walkers) - and is believed to represent an optimal 
foraging strategy [^. 

The dynamics of motile trajectories over shorter times, at the onset of dif¬ 
fusive dynamics, has received comparatively little attention but can inform on 
mechanistic aspects of the motile strategy. A prime example is the motility 
of E. coli bacteria, where three-dimensional tracking of motile cells revealed 
a “run-and-tumble” strategy of random motility, in which relatively straight 
uninterrupted paths (runs) are interspersed by rapid and random reorienta¬ 
tion events (tumbles) [^. The random walk of E. coli, therefore, is completely 
characterized by two random variables (run length and tumble angle) and two 
constant parameters (swimming speed and rotational diffusion coefficient). The 
effective diffusivity for its random translational motion can be expressed as a 
function of these parameters, and detailed studies over decades have yielded 
mechanistic models that link these key behavioral parameters to the underly¬ 
ing anatomy and physiology, including the sensory response to environmental 
stimuli mediated by a molecular signaling network 

Can a similar top-down approach be fruitfully applied to more complex 
organisms-for example, an animal controlled by a neural network? Animal 
behavior is both astonishing in its diversity and daunting in its complexity, 
given the inherently high-dimensional space of possible anatomical, physiologi¬ 
cal, and environmental configurations. It is therefore essential to identify appro¬ 
priate models and parameterizations to succinctly represent the complex space 
of behaviors — a non-trivial task that has traditionally relied on the insight 
and ingenuity of expert biologists. In this study, we ask if one can achieve a 
similar synthesis by an alternative, physically motivated approach. We seek 
a quantitative model with predictive power over behavioral statistics, and yet 
a parameterization that is simple enough to permit meaningful interpretations 
of phenotypes in a reduced space of variables. As an example, we focus on 
the motile behavior of nematodes (roundworms), which explore space using a 
combination of random and directed motility driven by undulatory propulsion. 
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The nematode C. elegans, long a model organism for the genetics of neural 
systems [^|^, is also ideally suited for quantitative studies of motile behavior. 
The worm’s behavioral repertoire [9 10 is commonly characterized in terms of 
forward motion occasionally interrupted by brief reversals [11[|12| , during which 
the undulatory body wave that drives its movement 13 switches direction. In 


addition, worms reorient with a combination of gradual curves in the trajectory 
(“weathervaning”) [14 15 and sharp changes in direction (“omega turns”) 11 
EH- These elementary behaviors are combined in exploring an environment 


14 


17 


14 


16 . Environmental cues such as chemical, mechanical, or thermal stimuli 


lead to a biasing of these behaviors, guiding the worm in favorable directions 
16[|18 . Finally, in practical terms, the worm’s small size (^1 mm in length), 
moderate propulsive speed (^100pms“^) and short generation time (^2 days) 
allow a considerable fraction of its behavioral repertoire to be efficiently sampled 
in the laboratory, on timescales of ^1 h within centimeter-scale agar arenas 
by video microscopy 


11 19 


20 . 


An influential example of such an analysis is the “pirouette” model proposed 
by Pierce-Shimomura and Lockery 16 which describes the worm’s exploratory 
behavior as long runs interrupted occasionally by bursts of reversals and omega 
turns that reorient the worm, in close analogy to the run-and-tumble model of 
bacterial random walks . Later work by lino et al. identified that worms also 
navigate by smoother modulations of their direction during long runs (“weath¬ 
ervaning”) [^, and recently, Calhoun et al. suggested that C. elegans may 
perform an infotaxis strategy in searching for food 21 . Importantly, while 


these previous studies have illuminated different modes of behavioral control, 
they were not designed to obtain a predictive model of the trajectory statistics 
and thus a succinct parameterization of C. elegans motility remains an impor¬ 
tant open problem. 

A quantitative parameterization capturing the repertoire of C. elegans' be¬ 
havioral phenotypes would facilitate data-driven investigations of behavioral 
strategies: for example, whether worms demonstrate distinct modes of motility 
(characterized by correlated changes in parameters) over time, across individuals 


in a population, or in response to changes environmental conditions 19 22 23 


A pioneering study on quantitative phenotypes has found that for morphology 
and gene expression, the distributions of naturally-occurring phenotypes are 
often constrained by evolutionary trade-offs between strategies represented by 
distinct parameter sets 24 . The motility of C. elegans and related nematodes 
provides a unique opportunity to apply an analogous approach to the study of 
behavioral strategies. C. elegans is a member of the Nematoda phylum, one of 



anatomy and ecological diversity makes nematode motility a compelling case for 
studies of behavioral phenotypes. Anatomical conservation suggests it might be 
possible to describe the behavior of diverse nematodes by a common model, 
and identifying the manner in which existing natural variation is distributed 
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across the parameter space of the model could reveal distinct motility strategies 
resulting from optimization under different environmental conditions. 

In this study, we develop a simple random walk model describing the trans¬ 
lational movements of a diverse collection of nematode species, freely-moving 
on a two-dimensional agar surface. In addition to providing a quantitative and 
predictive measure of trajectory dynamics, the parameters of our model define 
a space of possible behaviors. Variation within such a space can occur due to 
changes in individual behavior over time (reflecting temporal variation in the 
underlying sensorimotor physiology, or “mood”), differences in behavior among 
individuals (reflecting stable differences in physiology, or “personality”) and 
differences between strains and species (reflecting cumulative effects of natural 
selection). By quantitative analyses of such patterns of variation, we seek to 
identify simple, organizing principles underlying behavior. 


Results 


Nematodes Perform Random Walks Off-Food with a Broad 
Range of Diffusivities Across Strains 

In order to identify conserved and divergent aspects of motility strategies, we 
sampled motile behavior over a broad evolutionary range. We selected a phylo- 
genetically diverse collection of soil-dwelling, non-parasitic nematodes with an 
increased sampling density closer to the laboratory strain C. elegans (Figure]^ 
and Methods). To capture temporal and individual variation, we recorded the 
motility of up to 20 well-fed individuals per strain and each individual for 30 
minutes on a food-free agar plate. 

To enable measurements of behavioral statistics of multiple worms simulta¬ 
neously over long periods of time, we placed individual worms within a 10 x 10 
mm grid surrounded by a boundary of 1% sodium dodecyl sulfate (SDS), a 
detergent that C. elegans and most other nematodes avoided. While many C. 
elegans studies have used copper as a repellant boundary [^, we found that it 
did not strongly repel other nematodes (data not shown). 

We used automated image analysis techniques to identify worm outlines in 


each video frame (Figure S15). We then measured the centroid position {x{t)) 
and calculated the centroid velocity {v(t)). There was considerable variation in 
the spatial extent and degree of turning visible in the trajectories both within 
and across strains (Figure [^ . The boundary had subtle effects on be¬ 
havior, most noticeably a decrease in the extent of the trajectory on timescales 
beyond 100 s, but the statistics on shorter timescales were similar after excluding 
observations near the boundary (Figure [sTj). 

As previously seen in C. elegans 1^, the measured mean-squared displace¬ 
ment ([Aa;(r)]^) = {\x{t -|- r) — af(t)p) revealed a transition from ballistic to 
diffusive motion within a 100s timescale (Figure[I^, |S3[). Over short times, the 
worm’s mean-squared displacement scaled quadratically with the time lag r and 
speed s as (s)^t^ as expected for ballistic motion. Over longer times, the scaling 
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became linear with r as expected for random 2D diffusion (([Aa:(T)]^) = ADt). 
An effective diffusivity Dgff was estimated by fitting the observed mean squared 
displacements in the diffusive regime (see Methods). In agreement with the 
qualitative differences in the trajectories, the strains varied more than two or¬ 
ders of magnitude in this diffusivity (Figure [^, Table . 


The Random Walk of Nematodes Can Be Decomposed into 
Speed, Tnrning and Reversal Dynamics 

The broad range of observed diffusivities suggest that these diverse nematodes 
have evolved a variety of strategies for spatial exploration. To gain further 
insights into the manner in which such contrasting behaviors are implemented 
by each strain, we sought to further decompose the trajectory statistics with the 
aim of extracting a minimal model of the nematodes’ random walk. Whereas 
data from all ten measured strains were used in this development, for the sake 
of clarity in this and the following sections, we present our analysis and model 
development with a focus on data from three strains : CB4856 and PS312, which 
demonstrated two of the most extreme phenotypes, in addition to the canonical 
laboratory strain N2. (See SI for equivalent data for all strains.) 

The translational motion of the worm can be described by the time-varying 
centroid velocity v(t) which can in turn be decomposed into speed s{t) and 
direction of motion (hereafter referred to as its “bearing”) 

v{t) = = s{t) [cos(()(t),sin</)(t)] (1) 

To account for head-tail asymmetry in the worm’s anatomy, we additionally 
define the body orientation hereafter referred to simply as “orientation”) 

by the angle of the vector connecting the worm’s centroid to the head (Figure 
[1K)> identified using automated methods. The centroid bearing is related to 
this orientation of the worm by 




( 2 ) 


where the difference A'ijj(t) is a measure of the alignment of the direction of 
movement with the worm’s body orientation (hereafter referred to simply as 
“alignment”). We found for all strains that the distribution of was bi- 

modal with peaks at 0°and 180° (Figure}^,[S^). These match the forward and 


reverse states of motion described in C. elegans 11,12 


Each of the three components of the worm’s motility (speed, orientation, 
and alignment) varied considerably over time and in qualitatively different ways 
between strains (Figure]^). For example, the three strains shown in Figure]^ 
(N2, CB4856 and PS312) differed not only in their average speed, but also in 
the amplitude and timescale of fluctuations about the average speed. Similarly, 
a drift over time was obvious in the orientation of N2 but was smaller or absent 
in PS312, and the statistics of orientation fluctuations about the drifting mean 
also differed between strains. Finally, transitions between forward and reverse 
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runs were far more frequent in PS312 as compared to N2 and CB4856. Given 
the apparently random manner in which these motility components varied over 
time, we proceeded to analyze the dynamics of each of these three components 
as a stochastic process. 


Speed Dynamics 

Speed control has not been extensively studied in C. elegans, but it is known that 
worms move with a characteristic speed that is influenced by stimuli 17 . When 


intervals corresponding to transitions between forward and reverse runs were 
excluded from the time series, we found that the autocorrelation in speed fluc¬ 
tuations decayed exponentially over a few seconds (Figure [3|4 , [S5l4) , a timescale 
similar to the period of the propulsive body wave. These dynamics are natu¬ 
rally captured by an Ornstein-Uhlenbeck process [29| , which describes random 
fluctuations arising from a diffusive Wiener process (IFt) with magnitude \/2Ds 
that relax with timescale Tg back to an average value, /ig: 


ds{t) = Tg ^ [/Xg - s{t)] dt + \/2DsdWt 


( 3 ) 


Numerical integration of this equation closely reproduced the observed speed 
distributions during runs (Figure [^, |S5| 3). Transitions between runs were as¬ 
sociated with transient decreases in speed lasting a couple of seconds (Figure 
S6 4). The timescale of speed changes at these transitions was different from, 
and typically shorter than, the relaxation timescale of speed fluctuations during 


runs (Figure S6 3 


Diffusive Turning with Drift 

The orientation ijj(t) captures turning dynamics that are independent of abrupt 
changes in bearing due to reversals. To change orientation, C. elegans 
executes a combination of large, ventrally-biased “omega” turns 11 and 
gradual “weathervaning” [^. We found a weak bias in the worms’ orientation 
dynamics, resulting in a drift of magnitude (k^) of ^l°s“^ within each 100 s 
window (Figure , [S7]4) . 

After correcting for this drift, we found that the orientation correlation 
((cos['0(t) —'(/'(O)])) decayed nearly exponentially within a minute, consistent 
with previous studies in larger arenas [^. Thus, orientation fluctuations about 
the slow drift for each worm resembles a rotational diffusion process character¬ 
ized by a single, constant coefficient (Figure [^, [S7^). These observations 
lead to a simple model for the orientation dynamics that combines drift (approx¬ 
imated as a deterministic linear process over a 100 s window) with stochastic 
diffusion: 

d-ipit) = k^dt + ^/2D^dWt, (4) 


where dWt represents increments of a Wiener process 29 


We note that while the orientation drift in our data was well-described by 
a constant (k^) on the lOOs-timescale considered in our present analysis, both 
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the sign and magnitude of this drift did change over longer times (hundreds 
of seconds; Figur^S^. Such a slowly fluctuating drift could provide directional 
persistence over long times, as observed by Peliti et al. in C. elegans trajectories 
of durations > 1000 s We further note that the magnitude of in our 
data (~1 ° s“^) was similar to that of weathervaning excursions reported for C. 
elegans navigating in salt gradients 
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Forward and Reverse Runs 

The observation that motion switched abruptly between forward and reverse 
states that lasted for variable lengths of time suggested that reversals might 
be well described as a discrete random process. The simplest model for such 
a process is one in which transitions between the forward (A^/; = 0) and re¬ 
verse (A')/) = tt) states occur with a uniform probability in time, leading to 
exponentially distributed lifetimes for both forward and reverse runs: 


.^(^fwd run ^ ^/^fwd) 




>t)= exp{-t/Trev) 


(5) 

( 6 ) 


Such a two-state process (a ‘random telegraph’ process) yields an autocorrela¬ 
tion function that decays exponentially as = (cos(Ai/)(t-|-T) — A'0(t))) = 

CaiIj{oo) + (1 — Cav,(oo)) to a minimum value CAtp{oo) = ((rfwd — 

Indeed, we 


Trev)/(Trev + Tfwd))^ with a timescalc trt = + t ^ 2 ] 


31 


found that the Atjj autocorrelation function decayed nearly exponentially to a 
non-zero baseline (Figure]^, Figure S9 3). By fitting the autocorrelation func¬ 
tion with the two-state model, we extracted transition time constants that are 
in excellent agreement with the distribution of run durations beyond 1 s (Fig¬ 
ure]^). The shorter runs were not analyzed further as they appeared to occur 
mostly when the worm moved slowly, and therefore the bearing and velocity 
alignment was poorly defined (see example in Figure]^). 

In principle, the forward and reverse states could be characterized by dif¬ 
ferences in motility parameters other than these transition times, as forward 
and reverse motion are driven by distinct command interneurons in C. ele¬ 
gans 32p3 . However, we found that mean speeds were nearly identical between 
forward and reverse runs (Figure SIO). While we expect that this symmetry 


will be broken under some specific conditions, such as the escape response 34 


the strong speed correlation between the two states motivates the assumption, 
adopted in our model, that reversals change only the bearing (by 180°) and the 
propensity to reverse direction, represented in our model by the time constants 
'^fwd and Tj-ev- 
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A Model with Independent Speed, Turning and Reversals 
Captures the Ballistic-to-Diffusive Transition in Nematode 
Motility 


Given that the dynamics of the worm’s speed, turning and reversals could be 
described using simple stochastic processes, we asked whether combining them 
as independent components in a model of the worms’ random walk could suffi¬ 
ciently describe the observed motility statistics (Figure]^). We simulated tra¬ 
jectories of worms by numerically integrating equations (§- § for the speed, 
orientation, and reversal dynamics, respectively, which yields the worm’s veloc¬ 
ity dynamics through Eqs. Q and ([^, with Atp{t) equal to 0°during forward 
runs and 180°during reverse runs. Simulations of this model using parameters 
fit to individual worms produced trajectories that qualitatively resembled real 
trajectories and varied considerably in their spatial extent (Figure]^). 

Next, we quantitatively assessed the performance of the model in repro¬ 
ducing the statistics of the observed trajectories over the time scale of 100 s, 
within which all strains completed the transition from ballistic to diffusive mo¬ 
tion (Figure]^). As noted above (Figure]^) the collection of strains studied 
here demonstrated a broad range of effective diffusivities spanning « 3 orders of 
magnitude. We found that the model based on independent speed, turning and 
reversal dynamics closely reproduced not only the diffusivity of each strain but 
also the time evolution of the mean-squared displacement (([Aa;(T)]^)) across 
the ballistic-to-diffusive transition. (Figure]^, top). A closer inspection of the 
dynamics across this transition is possible by examining the velocity autocorre¬ 
lation function (C^ 7 (r)) the time integral of which determines the slope of the 
mean-squared displacement through {d/dt){[Ax{T)]‘^) = 2 dr'C^jlr'), a vari¬ 
ant of the Green-Kubo relation 35 ^ . The transition from ballistic to diffusive 
motion is characterized by the manner in which the normalized velocity auto¬ 
correlation Cij(T)/Cij{0) decays over the time lag r from unity (at t = 0) to zero 
(as T —7> oo). We found that these dynamics of C^{t) varied considerably across 
strains, not only in the overall ballistic-to-diffusive transition time tjj but also in 
the more detailed dynamics of the autocorrelation decay over time (Figure]^, 
middle). Salient features, such as the transition time, of the measured velocity 
autocorrelation functions G^j^obs were reproduced closely by the simulated ve¬ 
locity autocorrelation functions G^y^niodeh but there were also subtle deviations 
in the detailed dynamics for a number of strains. 

Given our model’s simplifying assumption that dynamics for s(t), ^nd 
A'0(t) are independent stochastic processes, we asked whether the remaining 
discrepancies between the simulated and measured velocity autocorrelation dy¬ 
namics could be explained by violations of this assumption of independence. 
As a model-independent assessment of the degree of non-independence, we first 
calculated the predicted velocity autocorrelation for the case that the dynam¬ 
ics of all three components are independent, G{ 74 ndep('r) = Gs(T)G,/,(T)GAi/>('r), 
where Gs(t), and are the autocorrelation functions of the mea¬ 

sured, rather than simulated, data for each of the components (see Supplemen¬ 
tal Info for derivation). We then compared the differences ~ G^j^ndep 




(blue curve in Figure |^, bottom) and Cij^obs ~ C'ff,model (red curve in Figure 
1^, bottom). Indeed, there were subtle differences both on shorter (~ls) and 
longer timescales (~10s). However, where differences did occur for the simu¬ 
lated model, the magnitudes of differences were very similar to, or less than, 
those for the product of component-wise measured correlation functions (ie., 
C'u.obs ~ C'u,model ^ C'-D’.obs ~ C'lT.indep)- Wc therefore conclude that the relatively 
subtle differences between the data and model arise primarily in instances where 
this assumption of independence between the three motility components (speed, 
turning, and reversals) breaks down. 


Low-Dimensionality in Behavioral Variation across Species 


The results presented in the previous sections demonstrate that a random-walk 
model with a total of seven parameters describing independent speed, turning 
and reversal dynamics, provides a good approximation of the worms’ motile be¬ 
havior. The model parameters thus define a seven-dimensional space of motility 
phenotypes in which behavioral variation across strains and species can be ex¬ 
amined. If components of behavior were physiologically regulated or evolution- 
arily selected for in a coordinated manner, we would expect to find correlated 
patterns in the variation of these traits. 

We fit our model to the trajectory statistics of each individual worm and 
built a phenotypic matrix of 139 worms x 7 behavioral parameters (for a sum¬ 
mary of all parameters for each strain, see Tables S2|S4 ). To characterize the 
pattern of variation within the phenotype space, we computed the covariance 
matrix for the seven parameters (Figure [sjA). The parameters with the largest 
variance (diagonal elements in Fig. 5A) were the forward and reverse state life¬ 
times (r/ujd, Trev), followed by those describing speed dynamics (ts, Ds). More 
interestingly, there was extensive covariation among the model parameters, not 
only within the parameters of each motility component (speed, orientation, re¬ 
versals) but also between those of different components. 

We looked for the dominant patterns in the correlations using principal com¬ 
ponent analysis (Figure]^), a technique which identifies linear combinations of 
parameters that represent orthogonal axes (modes) capturing the most variance 
in the data [37| . The distribution of variance captured by these modes uncovered 
a single dominant mode of correlated variation (Figure]^, left). This principal 
mode, capturing about 60% of the total variation, described significant correla¬ 
tions among all the parameters except for Dg (Figure right. Table S5). We 
did not attempt to interpret higher modes since, individually, none significantly 
exceeded the captured variance under a randomization test (see Methods, and 
Figure]^, left). A non-linear technique for describing the behavioral manifold 
did not reveal any distinct additional structure (Figure [Sll| ). 

We used numerical simulations to determine the effects on motile behavior of 
varying parameters along the principal mode. The measured trajectory pheno¬ 
types projected onto this mode in the range {—2, 2} centered about the average 
phenotype located at the origin, and we performed simulations for parameter 
sets evenly sampled along this range. As the projection on the principal mode 
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increased, the simulated trajectories became more expansive and the effective 
diffusivity increased by more than two orders of magnitude (Figure]^). This 
largely reproduced the observed variation in the measured diffusivities Deff as 
a function of the projection along the first mode. 

The significant loadings on nearly all of the model parameters indicated that 
this dominant mode of behavioral variation reflects global changes in motil¬ 
ity. The manner in which the motility parameters change along this mode 
are reminiscent of conventional descriptions of roaming and dwelling behav- 
and other organisms 


lor m worms 


19 


38 39 : Roaming worms move at high 


speed with infrequent reversals, while dwelling worms move at low speed and 
reverse frequently 19 . Interestingly, the parameters with significant loadings 


on this mode included also those pertaining to gradual changes in orientation 
(fc^, Il^),and the sign of these loadings indicate that roaming worms tend to 
generate trajectories with higher curvature than do dwelling worms. Changes 
in path curvature have not been, to our knowledge, implicated in previous stud¬ 
ies of roaming/dwelling in nematodes and hence could represent a previously 
unrecognized difference between these behaviors. 


Specialized and Diversified Behavioral Strategies Across 
Strains 

The principal behavioral mode discussed in the preceding section was identi¬ 
fied by analyzing variation across all individual worms measured in this study. 
However, these worms came from diverse strains and species that differ in their 


average behavior (see Tables SI - S4). How does the variability among individ¬ 


uals of a given strain compare to differences between the average phenotypes of 
strains/species? On the one hand, each strain might be highly “specialized”, 
with relatively small variation within strains as compared to that across strains. 
On the other hand, strains might implement “diversified” strategies in which 
genetically identical worms vary strongly in their behavior. To address these 
two possibilities, we analyzed the distribution of individual phenotypes within 
each strain, as well as that of the set of averaged species phenotypes. 

For each measured individual, we computed the projection of its motility 
parameter set along the principal behavioral mode and estimated strain-specific 
distributions of this reduced (one-dimensional) phenotype (Figures |S12[ Ta¬ 
ble S6). While the precise form of these distributions could be relevant for 


evolutionary fitness, we focused our analysis only on their first two moments, 
characterized by the mean and standard deviation, given the finite number of 
individuals (< 20) per strain. Further, we computed the principal-mode projec¬ 
tion of the average phenotype of each species to define an interspecies phenotype 
distribution (Figure]^. 

Strains varied considerably in both the mean and standard deviation of 


their phenotypic distributions (Figure S13). Remarkably, the standard devia¬ 
tion for individual variation within each strain was comparable in magnitude 
to that for the set of average phenotypes across species. Some strains were 
specialized towards roaming or dwelling behavior, such as CB4856 and PS312, 


10 














respectively, with a strong bias in their behavior and comparatively low indi¬ 
vidual variability. Others, such as N2 and QX1211, appeared more diversified 
with an intermediate average phenotype and higher individual variability. These 
considerable differences in phenotype distributions across strains reveal the evo¬ 
lutionary flexibility of population-level heterogeneity in nematodes, and suggest 
a possible bet-hedging mechanism for achieving optimal fitness in variable en¬ 
vironments 


40,41 


In assessing such variability of phenotypes, it is essential to ask how uncer¬ 
tainty in the determined parameters (obtained from model fits) contribute to 
the observed variability in phenotypes. We therefore computed the contribu¬ 
tion of uncertainties in the individual-specific fits of all parameters, which sets a 
lower bound on such uncertainties, by way of error propagation, (see Methods) 
This analysis revealed that parameter uncertainty contributed less than 10% 
on average to the variation in the reduced phenotype across the entire dataset, 
and was also many-fold lower than the individual variation within each strain 
(Figures!^ |S12[ ). Furthermore, strains that were closer on the phylogenetic tree 
tended to have more similar phenotypic distributions (Figure [S14[ ). These find¬ 
ings support the view that the phenotypic variation estimated in the current 
analysis largely represented true differences in behavior. 


Discussion 


We have presented a first systematic analysis of motile behavior across a broad 
range of strains and species of the nematode phylum, ranging from the lab 
strain C. elegans N2 to Plectus sjh2 at the base of the chromadorean nematode 
lineage. Despite the large evolutionary distance (in terms of ribosomal RNA, 
the divergence of species studied here are comparable to that of humans from 


fish 42 ), we found that a behavioral model described by seven parameters 
could account for much of the diversity of the worms’ translational movement 
across the ~ lOOs-timescale spanning the ballistic-to-diffusive transition. This 
simple model provides a basis for future studies aiming to capture more detailed 
aspects of nematode behavior, or to connect sensory modulation of behavior to 
the underlying physiology. More generally, this study shows that quantitative 
comparisons of behavioral dynamics across species can provide fundamental 
insights regarding the design of behavioral strategies. We summarize below our 
main findings and discuss their implications, highlighting along the way open 
problems that represent promising directions for future work. 


The Minimal Model: What Does It Capture, and What 
Does It Miss? 

In this study, we focused on a high-level output of behavior — translational 
and orientational trajectory dynamics — and sought to build the simplest pos¬ 
sible quantitative model that could capture the observed behavioral statistics. 
We found that a model with only three independent components — (I) speed 
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fluctuations that relax to a set point on a timescale of a few seconds, (2) ori¬ 
entation fluctuations with drift, and (3) stochastic switching between forward 
and reverse states of motion — describes well, overall, the trajectory statistics 
of all tested nematode species across the ballistic-to-diffusive transition in the 
absence of stimulus (Figure]^. 

Notably, we have not included explicit representations of some reorienta¬ 
tion mechanisms that have been studied in the past. One such mechanism is 
the “pirouette”, which has been described as a mechanism for abrupt reorien¬ 
tations in C. elegans. Pirouettes were initially identified from biexponentially 
distributed “runs” separated by either changes in body wave direction or sharp 
In our data, we find that the timing of the initiation and termination 
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turns 

of reversals, which would both count as runs in the pirouette description, follow 
exponential distributions with similar time constants as previously reported for 
the pirouette run distribution. A second reorientation behavior not explicitly 
represented in our model is the “omega” turn, a sharp change in direction in 
which the worm’s head touches its body [^. While omega turns might in¬ 
deed be mechanistically distinct from gradual turns, their contribution to the 
trajectory statistics within our data set were likely modest, as we were able to 
adequately describe turning within our data set by a continuous diffusion-drift 
process (Figure [^,D). It is possible, however, that explicit representations of 
pirouettes and/or omega turns would be important in other experimental sce¬ 
narios, e.g. those that include navigation in the presence of gradient stimuli. 

In its current form, our simple model does not account for possible cor¬ 
relations between the dynamics of the three motility components (speed, ori¬ 
entation, and reversals). Comparisons of simulated versus measured trajecto¬ 
ries demonstrated that the effects of such correlations on the motility statis¬ 
tics are small but detectable (Figure]^). The differences were most significant 
for the velocity-autocorrelation dynamics on a ~10s timescale, and were simi¬ 
lar to those for model-free predictions obtained by combining component-wise 
correlation functions under the assumption of independence. Discrepancies on 
this intermediate timescale occurred most often in fast-moving strains (CB4856, 
JU775, sjh2) that frequently approached the repellent boundary. Therefore, we 
suspect that the discrepancy arises from a behavior such as the escape response 
that has been described in C. elegans, in which worms respond to strongly aver¬ 
sive stimuli by reversing and then making a sharp omega turn over the course of 
several seconds [34]. Such a stereotyped sequence would be expected to intro¬ 


duce temporal correlations between speed changes, turning, and reversals that 
would be missed in our model with independent speed, orientation and reversal 
dynamics. 

While here we have focused on the transition to diffusive motion, some recent 
experiments suggest that C. elegans might engage in superdiffusive behavior on 
timescales longer than 100 s 15 23 . Superdiffusive behavior could arise from 
nonstationarities in motile behavior on these timescales, such as behavioral state 

Another mechanism for 
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transitions occurring on a timescale of minutes 
superdiffusion is directed motility 
chemical or thermal gradients. In such environments 


in response to external stimuli such as 
nematodes are known 
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to use at least two distinct mechanisms for navigation 14 16 and the model 


here could be extended by studying the dependence of motility parameters on 
environmental statistics. 

Another promising approach for building a more complete behavioral model 
is to incorporate information about the body shape 28,44 . Centroid behavior 


can confound or hide distinct dynamical patterns in postural motion; for exam¬ 
ple, speed can be controlled both by the shape and frequency of the locomotory 
body wave. Indeed, recent work by Brown et al. showed that a rich repertoire 
of dynamics can be identified as temporal “motifs” in the postural time series 


of C. elegans and used to classify mutants with high discriminatory power 45 


We have found that all of the species tested here can also be described with a 
common set of postural modes (not shown), suggesting future directions on the 
evolutionary space of postural dynamics. 


The Roaming/Dwelling Behavioral Mode: Variability and 
its Physiological Basis 


While we found that a single behavioral model could be used to characterize 
nematode motility across the chromadorean lineage, the parameters of the model 
varied extensively from strain to strain. Quantitatively, more than half of the 
variation corresponded to a correlated change in the parameters underlying 
the timing of forward and reverse runs and the dynamics controlling speed and 
turning (Figure]^). This pattern of parameter modulation drove a change from 
low speed, short, straight runs to high speed, long, winding runs and is similar 
to the canonical definition of roaming and dwelling in C. elegans 
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Roaming and dwelling are thought to represent fundamental foraging strate¬ 
gies reflecting the trade-off between global exploration and local exploitation, 
respectively, of environmental resources 46 . Recent work has suggested that 


such archetypal strategies can be recovered by quantitatively analyzing the ge¬ 
ometry of phenotypic distributions in parameter space 24 . The motility pheno¬ 


types we found in the present study were biased along one principal dimension, 
with the extremes corresponding to roaming and dwelling behaviors. This ob¬ 
servation compels us to suggest that an exploration-exploitation trade-off is the 
primary driver of phenotypic diversification in the motility of chromadorean ne¬ 
matodes in the absence of stimuli. Interestingly, a recent study on the motility 


of a very different class of organisms (ciliates) yielded a similar conclusion 39 
across two species and different environments, the diversity of motility pheno¬ 
types was found to be distributed principally along an axis corresponding to 
roaming and dwelling phenotypes. The emergence of roaming/dwelling as the 
principal mode of variation in such disparate species underscores the idea that 
the exploration-exploitation trade-off is a fundamental constraint on biological 
motility strategies. 

A surprising finding in our study was that the extent of behavioral variabil¬ 
ity across individuals within a strain was comparable to that for variation of 
phenotypes across species (Figure [6 S12). In slowly changing environments, 
the most evolutionarily successful species are those that consistently perform 
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well in that environment. This can be achieved by evolving a specialized, high 
fitness phenotype that varies little among individuals. However, increased phe¬ 
notypic variability among individuals can improve fitness in more variable envi¬ 
ronments if some individuals perform much better in each condition-a so-called 
“bet-hedging” strategy 40 41 . The large variability we observed among indi¬ 


vidual phenotypes within each strain might reflect such a bet-hedging strategy 
in nematode exploratory behavior. 

The observation that individual variation is comparable to inter-species 
raises the intriguing possibility that there exist conserved molecular and/or 
physiological pathways driving diversification of spatial exploration strategies. 
Analogous variation in exploratory behavior was also detected in a recent analy¬ 
sis of nonstationarity in the behavior of wild-type and mutant C. elegans under 
various nutritional conditions [^. Physiologically, protein kinase G (PKG) sig¬ 
naling and DAF-7 (TGF-/3) signaling from the AST neuron are thought to be 
major mechanisms controlling roaming and dwelling in C. elegans [19[|22| . In- 
triguingly, PKG signaling is also involved in controlling foraging in Drosophila 
and other insects as well as many aspects of mammalian behavior 47,48 . Re¬ 


cently, Flavell et al. also elucidated a neuromodulatory pathway involving sero¬ 
tonin and the neuropeptide pigment dispersing factor (PDF) controlling the 
initiation and duration of roaming and dwelling states [43| . 

Perturbations to the molecular parameters of such pathways underlying 
global behavioral changes might provide a mechanism for the observed cor¬ 
related variations at the individual, intra-, and inter-species levels. Our simple 
model provides a basis for future investigations to uncover conserved mecha¬ 
nisms that generate behavioral variability, by dehning a succinct parameteriza¬ 
tion of behavior that can be combined with genetic and physiological methods. 
The identification of such conserved pathways affecting many phenotypic pa¬ 
rameters is of fundamental interest also from an evolutionary perspective, as 
they have been proposed to bias the outcome of random mutations towards 


favorable evolutionary outcomes 49 50 


Materials and Methods 


Selection of Strains 

The nematode phylum is classically divided into three major branches—chromadorea. 


enoplea, and dorylaimia—that are broken into a total of five major B-clades 51 


and twelve minor H-clades 52 . The chromadorean lineage is the largest, span¬ 
ning B-clades III-V and H-clades 3-12 [^[^ . C. elegans is located in clade V9 
(the rhabditids), one of the most diverse clades 53 . In addition to the lab strain 


N2, we selected three of the most genetically distinct wild isolates of C. elegans 
(GB4856, JU775, and QX1211) to sample intraspecies variation 54 . From H- 


clade 9 in order of increasing evolutionary distance, we selected Caenorhabditis 
briggsae JU757, Diploscapter sp. PS2045, Rhabditis myriophila DF5020, and 
Pristionchus pacificus PS312. The next closest major group, B-clade IV, con- 
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tains H-clades 10-12. H-clade 12 contains the plant parasitic tylenchs and was 
thus not included in this study. H-clades 10 and 11 contain many bacterial feed¬ 
ers, of which we selected Panagrolaimus sp. PS1159. Finally, from the basal 
chromadorea, we obtained Plectus sp. sjh2, a member of H-clade 6 . 

C. elegans N2, CB4856 and JU775 were provided by the Caenorhabditis 
Genetics Center, which is funded by NIH Office of Research Infrastructure Pro¬ 
grams (P40 OD010440). C. elegans QX1211 was kindly provided by Erik An¬ 
dersen (Northwestern Univ.). Plectus sp. sjh2 was isolated from a soil sample 
using morphological criteria by Casper Quist and Hans Helder (Wageningen 
Univ.). SJH then isolated a single species by starting cultures with a single 
worm. The remaining strains were used in previous studies by LA [55] . 


Cultivation of Worms 


Worms were grown on NGM-SR plates (3g NaCl, 24 g agar, 2.5 g peptone, 
ImL 5mgmL“^ cholesterol in EtOH in 975mL water, with 1 mL iM CaCU, 
ImL 1 m MgS 04 , 25mL 1 m K 2 PO 4 pH 6 , ImL 200mgmL“^ streptomycin in 
water, and 0.23 g 5mL 40mgmL“^ nystatin in DMSO (added after autoclaving) 
seeded with E. coli HBlOl, as previously described [^. E. coli HBlOl was 
first cultured in M9 minimal media (3g KH 2 PO 4 , 6 g Na 2 HP 04 , 5g NaCl, ImL 
1 m MgS 04 in IL water) supplemented with 10% Luria broth and lOmgmL”^ 
streptomycin [^. Plates were incubated with a light circle of HBlOl culture 
for a day at 37 °C and then stored at 4°C. For Plectus sp. sjh2, low salt plates 
(2% agar supplemented with 5mgL“^ of cholesterol from a 5mgmL“^ EtOH 
solution) were used as previously described 58 . On NGM-SR plates, these 


worms became shriveled and died. As the plates did not have nutrients for the 
bacteria to grow, HBlOl was grown to high density in Luria broth overnight at 
37 °C, washed 3X in water, resuspended at lOX concentration, and applied to 
the plates. 

Nematodes were cultured by either transferring a few worms by worm pick or 
a chunk of agar to a new plate after the worms reached adulthood. The plates 
were then incubated at 20 °C. The growth rate varied considerably among 
strains, with Plectus sp. sjh2 taking nearly two weeks to reach adulthood. We 
avoided starving the worms at any point during their cultivation, especially in 
the period before behavioral experiments were performed, as this can induce 
transgenerational phenotypic changes 59,60 , and we have observed transient 


effects on motility lasting at least a couple of generations (data not shown). 


Imaging 

The imaging experiments were done on 3.5 cm plates containing the same me¬ 
dia used for cultivation. A 2x2 10 mm repellant grid was made by etching the 
plate with a tool dipped in 1% sodium dodecyl sufate. Four young adult, well- 
fed nematodes were transferred individually by worm pick to a 10 pL drop of 
M9 (water for Plectus sp. sjh2) to remove bacteria stuck to the worms. The 
worms were then transferred by pipette in a minimal amount of buffer to the 
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imaging plate, and excess buffer was removed as much as possible. The plate 
was imaged 10-20 minutes after picking the worms, minimizing most transient 
behaviors. The plate was placed on a custom imaging rig in an inverted, uncov¬ 
ered configuration with illumination by a Schott MEBL-CR50 red LED plate. 
The behavior was recorded for 30 minutes using a Point Grey Grasshopper Ex¬ 
press GX-EW-60S6M-G camera equipped with an Edmund Optics NT54-691 
lens (set to a magnification of 0.5X) at a resolution of 2736x2192 (12.5pm/px) 
at 11.5 frames/s using a custom National Instruments LabView acquisition pro¬ 
gram. The video was subsequently compressed using the open-source XVid 
MPEG-4 compression algorithm using maximal quality settings. 


Tracking and Image Analysis 


The behavioral videos were analyzed using a custom automated analysis pro¬ 
gram in MathWorks Matlab. The average background was calculated from 50 
frames evenly sampled across the entire video. The background was then sub¬ 
tracted from each frame and a global threshold was applied. The thresholded 
image was cleaned by applying a series of morphological operations: Incomplete 
thresholding of the worm was smoothed by applying morphological closing with 
a disk with a similar radius as the worm. Any remaining holes were filled in 
using a hole-filling algorithm. Small holes or ones with a low perimeter to area 
ratio were excluded as they sometimes fill in worms undergoing an omega turn, 
as described in 61 . Finally, regions in which the worm was just barely touching 


itself were split by sequentially applying open, diagonal fill, and majority mor¬ 
phological operations. The worm was then identified as the largest connected 
component with an area within 2-fold of the expected value. The centroid was 
used for tracking. In addition, the image skeleton was calculated. Sample im¬ 
ages from each of the processing steps are shown in Figure [S15| 

The head of the worm was automatically identified from the statistics of the 
worm’s posture. The skeleton in each image frame was fixed by minimizing the 
total distance between the skeleton points in two possible orientations with the 
last identified skeleton. The trajectory was then divided into segments of at 
most 500 frames containing no more than ten consecutive missing skeletons and 
with at most an average distance of 10 pm per skeleton point. Segments with 
less than 150 frames were discarded. For each segment, the head was assigned 
by comparing the summed variance in body angles of the front and back 10% 
of the skeleton. The end with the greater variance was selected. 

The velocity v(t) was calculated from the centroid position v(t) using the 
derivative of a cubic polynomial fit to a sliding 1 s window. The direct esti¬ 
mation of the velocity using a symmetrized derivative had a large (5-correlated 
component that interfered with later analysis. The use of the cubic polynomial 
did not noticeably distort the correlation functions (Figure [Sl^. 
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Calculation of Behavioral Statistics 


The worm’s behavior fluctuated or sometimes drifted over long times (Figure 


S4), but the average statistics over 100 s windows were approximately stationary. 


In order to focus on dynamics within the 100 s timescale, the mean-squared 
displacement and all auto- and cross-correlation functions were calculated for 
100 s windows and then averaged. This reduced the influence of longer timescale 
fluctuations in the speed and reversal rate. For all calculations, observations 
near the boundaries and pairs of points between which the worm approached 
the boundary were excluded. 

The uncertainty of the motility parameters and behavioral mode projection 
of individual worms were estimated from the uncertainty of the curve fitting 
(described in the subsequent sections). The uncertainty of the parameters of 
each curve fit was calculated from the estimate of the Jacobian matrix returned 
by the curve fitting routine using the nlparci tool in MathWorks Matlab. The 
standard deviation a of each curve fitting parameter was then estimated from 
the 95% confidence interval Cl as ct « Cl /4. Finally, these standard deviations 
were used with standard error propagation formulas to estimate the uncertainty 
of the motility parameters and behavioral mode. 


Calculation of Effective Diffusivity, /Jeff 

To estimate the effective diffusivity Ueffj we fit the mean-squared displacement 
([Aa;(r)]^) over the diffusive regime. For this purpose, we defined the diffusive 
regime as the time-lag interval after which the normalized velocity autocorrela¬ 
tion Cj;{t)/Cj;{Q) decayed to 0.1. We note that in some cases (especially for fast- 
moving strains such as CB4856, JU775 and sjh2) the fit to ([Aa;(T)]^) = dUeffT 
in this regime was poor due to boundary effects arising from the finite size of the 
behavioral arena. For these strains, Deff should be regarded as a lower bound 
for the true diffusivity. 


Reversal Analysis 

The reversal state was assigned as described in the main text using the A'0(t) 
variable. Runs were identified from contiguous segments of forward or reverse 
motion. Frames with missing postural data were marked as ambiguous ends of 
runs. Assuming a stochastic telegraph process that generates states A')/) = 0 
(forward) and /S.ip = tt (reverse) with probabilities 1 — frev and frev, respectively, 
the autocorrelation at long time lags is Ca,Ij{t —>■ oo) = (1 — 2frevY- For the 
proposed telegraph process, each state has an exponentially distributed lifetime 

^r6v 

(Tfwd, TVev) and therefore frev = -■ The expected correlation timescale 

Trev + "rfwd 

for the mixture of the two states is g(rrev, Tfwd) = (Rwd + Rev) ■ The Atp 
autocorrelation function was therefore fit to 


Caip{t) 


[1 - /(TrevjTfwd)] SXp 


T 

'^fwd) 


F /('^rev,'n^wd) 


(7) 
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where /(rrev,Tfwd) = The fraction of time spent reversing is: 

^fwd “t ^rev 

/rev = 0.5—-^/(Trev, T-fwd)/4, where /rev S [0, 0.5]. The transition time constants 

^(S^revj "^fwd) , ^/(Tev; Twd) 

and rtwd = 


are then Trev = 


1 -/. 


/r. 


Speed Analysis 

Transitions between forward and reverse runs were ignored in the speed data 
by excluding data points within 1 s of the start or end of a run. The speed 
set point fj,s was fit by taking the mean. The remaining parameters of the 
speed dynamics ([^ were fit by its analytical autocorrelation function: Cs(r) = 
DsTs exp {-t/ts). 


Orientation Analysis 

Drift in the body bearing (V'(t)) was identified by fitting the linear trend in 
the unwrapped bearing over 100s windows, ip(t) = + ipQ. The average 

absolute value of the drift rate, (|fc^|), across the ensemble of windows was used 
to characterize the trajectory as a whole. The autocorrelation function was then 
fit to an exponential as described in the text to obtain D^. 

Simulations 

Reversals, orientation, and speed dynamics were all simulated independently 
using the model described. Forward and reverse run durations were chosen 
according to equations © and © by drawing exponential random numbers with 
mean value rf^d or Trev During reverse runs, A'lp was set to tt. The orientation 

g and speed ([^ dynamics were simulated using the Euler-Maruyama method 
with a time step that matched the frame rate. To prevent negative speeds, 
a reflective boundary condition was imposed by taking the absolute value of 
the speed at each simulation step. The velocity was then calculated from the 
decomposition in Q and trapezoidally integrated to give the centroid position 
x{t). 

Behavioral Mode Analysis 

The model parameters were ht to each trajectory to give a phenotypic matrix 
T. The phenotypic matrix was centered by subtracting the mean phenotype, 
T = T — (T)indiv.- The covariance matrix was then calculated, Cx = covT, 
and decomposed into eigenvalues A and eigenvectors (behavioral modes) b, 
Cxb = Ab. To reduce any bias coming from a single trajectory or strain, this 
calculation was bootstrapped lOOOX by simultaneously replacing both trajecto¬ 
ries within each strain and a single strain with replacement. The significance 
of each mode was assessed by comparison with the eigenvalues of samples of 
a random matrix model in which correlations between model parameters were 
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removed. Correlations between model parameters were removed by indepen¬ 
dently shuffling observations within each of the seven columns of the original 
139 X 7 matrix. The projections of each trajectory on these behavioral modes 
were calculated by P = Tb. 


Statistics 


Unless otherwise indicated, errorbars represent a 95% confidence interval es¬ 
timated from 1000 bootstrap samples. All probability distributions were em¬ 
pirically estimated using kernel density methods in MathWorks Matlab with a 
bandwidth automatically selected using Silverman’s rule of thumb 
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Figure 1: Nematodes perform random walks off-food with an effective diffu- 
sivity that varies across strains. (A) A phylogenetic tree with the strains used 
in this study is shown. The bold numbers are the major clades of Nematoda. 
The gray box indicates genetically distinct wild isolates of C. elegans. A picture 
of a worm and an individual 30 minute trajectory are shown to the right. In 
addition to inter- and intra-species variation, up to 20 individual worms were 
tracked for up to 30 minutes, providing additional information about individual 
and temporal variation. (B) The average mean-squared displacement across N2 
individuals is shown in black. For comparison, the expected ballistic (blue, cal¬ 
culated from the average speed) and diffusive (red, fitted to the diffusive regime; 
see Methods) motility are plotted. Consistent with previous reports, the motil¬ 
ity transitions from ballistic to diffusive over 100 s. The shaded regions indicate 
95% confidence intervals. (C) Effective diffusivity (mean and 95% confidence 
intervals) for each strain, calculated from fits of the mean-squared displacement 
as in B. 
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Figure 2: The random walk of nematodes is composed of speed, turning, and 
reversal dynamics. (A) The motility of the worm is described by a time-varying 
centroid position x{t) and physical orientation '0(t). From this, the speed s{t) 
and a measure of the alignment of the velocity with the orientation Atplt) are 
calculated. (B) One minute of speed, orientation, and velocity alignment data 
from the 30 minute experiment is shown for individuals from three strains. 
(C) The probability distribution of the velocity alignment reveals a bimodal 
structure corresponding to forward and reverse motion states. Shaded areas 
represent 95% confidence intervals. 
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Figure 3: Statistical characterization of the dynamics of the motility compo¬ 
nents. (A-B) The speed fluctuations excluding reversal state transitions were 
characterized by an Ornstein-Uhlenbeck process. (A) The autocorrelation of the 
speed indicated that fluctuations decayed exponentially over a few seconds. (B) 
The model reproduced the observed speed distributions. (C-D) The orientation 
dynamics were modeled with drift and diffusion. (C) The mean absolute change 
in orientation drifted linearly with time. (D) The orientation autocorrelation 
function decayed exponentially over 100 s as expected for diffusion. (E) The 
velocity alignment autocorrelation decayed over 10 s to an intermediate value. 
Similar behavior is expected for a simple two-state Poisson process model in 
which worms transition between forward and reverse states with exponential 
waiting times. (F) The cumulative distribution of the forward and reverse run 
durations for an individual N2 worm are plotted on a log scale. There were a 
large fraction of short events on a less than 1 s timescale followed by an expo¬ 
nential regime. The time constants estimated from the fit of the autocorrelation 
function in E were used to plot the expected exponential distribution after fit¬ 
ting the observed fraction of short events (red). In each plot, the ensemble 
average of individuals from the strains are shown with solid lines. Shaded areas 
represent 95% confidence intervals. The model fits are shown in dashed lines. 
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Figure 4: An independent model of speed, turning, and reversals quantita¬ 
tively captures nematode motility. (A) Summary of the model. (B) Simulated 
trajectories for the three example strains. (C) Statistical comparison of the 
data (black) and simulations (red), ensemble averaged across individuals for 
each strain. The mean-squared displacement (top) was closely reproduced in 
all cases. The normalized velocity autocorrelation, C's(r)/C'^T(0), (middle) was 
less well captured. However, the errors in the velocity autocorrelation (bottom) 
from the simulations (red) arise from assuming the dynamics of the speed, ori¬ 
entation, and velocity alignment are independent (blue). Shaded areas indicate 
the 95% confidence interval. 
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Figure 5: Motility parameters co-vary along an axis controlling exploratory 
behavior. (A) Covariance matrix of the behavioral parameters across the whole 
dataset. The diagonal indicates the variance for each parameter, while off- 
diagonal elements indicate covariation among parameters. (B) (left) Significant 
patterns (eigenvectors) in the covariance matrix were identified by comparing 
the fraction of variance captured by each with the amount expected for a un¬ 
correlated dataset (red line), (right) The loadings on the top eigenvector (be¬ 
havioral mode). (C) The effective diffusivity (left, circles) and a 30 minute 
trajectory (right, colors match points on graph) from simulations in which the 
average phenotype was varied along mode 1. The projections and effective diffu¬ 
sivity of the measured trajectories are shown as black points, and the average of 
each strain is shown as a gray square. Error bars and shaded regions represent 
bootstrapped 95% confidence intervals. 
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Figure 6: Specialized and diversified behavioral strategies across strains. Phe¬ 
notypic distributions showing the variability among individual trajectories sam¬ 
pled from a single strain (N2, CB4856, and PS312) and (bottom) the variability 
in the mean phenotype of each species are shown. Red dashes indicate individ¬ 
ual observations, the distributions shown are kernel density estimates. The 95% 
confidence interval for a single, randomly chosen data point is shaded in gray 
for reference. 
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Supplemental Information 

Derivation of the Velocity Autocorrelation Function Under 
the Assumption of Independence 

The velocity autocorrelation function can be written in terms of the motility 
components, 

Cy{T) = (u(0) • u(t)) 

= (s(0) [cos [t/’(0) + A'0(O)], sin [^/>(0) + A^/;(0)]] 

s(r) [cos [ipir) + AiP{t)] , sin [-ipir) + A'(/'(t)]]) 

The expected value of the product of independent random variables is the prod¬ 
uct of the expected value of each variable, i.e. {xy) = {x){y). Therefore we can 
factor out Cg = (s(O)s(t)), leaving the vector product with tp and The 

expanded vector product is: 

Cv{t) = Cs(t)x (cos [V'(O) -I- Ai/)(0)] cos [V'(t) -f A'!/)(t)] -fsin [^’(0) -I- A^(0)] sin + A'0(t)]) 

( 8 ) 

The trigonometric functions on '0(t) + A'0(t) can be rewritten as products of 
trigonometric functions of the terms: 

cos [V'(t) -I- A^(t)] = cos ^(t) cos Aip{t) — smtp{t) sin A'tjj{t) 

sin [ijj{t) + Atp{t)] = sm'ijj{t)cosAtl}{t) + costp{t)sinA'ip{t) 

However, since Aip{t) = {0,7r}, sinAip^t) = 0: 

cos[4’{t) + Ail;{t)] = cos tp{t) cos A'tjj{t) 
sin ['0(t)-I-AV'(t)] = sin'0(t) cos A^(t) 

Substituting into ([^, 


Cv{t) = Cs(t)x (cos V'(O) cos 4>{t) cos A'ip(O) cos A-iP(t)+ sin'0(0) sin'0(r) cos A'0(O) cos A0(t)) 

( 9 ) 

We can now factor out C^{t) = (cos [0(t) — 0(0)]) = (cos0(0) cos 0(t) -f 
sin0(O)0(T)) to get: 

Cv{t) = C's(t)C',/,(t)(cos A0(O) cos A0 (t)) 

Finally, we substitute (again dropping sin A0(t) terms): 

CaiIj{x) = (cos [A0(O) — A0(r)]) = (cos A0(O) cos A0(t)) 


to get: 


C^T,indep(r) = Cs{t)C^{t)Ca^{t) 
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Figure SI: The boundary does not substantially effect the ballistic to diffu¬ 
sive transition. Comparison of the statistical behavior of C. elegans N2 in the 
experiments presented here with a 1% SDS boundary (black) and a previously 
reported dataset 28 (red). The mean squared displacement and normalized 
velocity autocorrelation, Cij{T)/Ci;{0), are shown, as in Figure]^ Shaded areas 
indicate the 95% confidence interval. 
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Figure S2: Overview of the dataset. Trajectories of all worms, tracked for 30 
minutes, included in the study. The axes correspond to the 10 mm by 10 mm 
boundary. Points excluded from the analysis because they were close to the 
boundary are highlighted in blue. 



Figure S3: Ballistic to diffusive transition for all strains. The average mean- 
squared displacement (MSD), calculated across individual trajectories, is shown 
in black for each strain. The expected ballistic and diffusive MSD curves, cal¬ 
culated as in Figure [^, are shown in blue and red, respectively. Shaded areas 
indicate the 95% confidence interval. 
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Figure S4: Average changes in behavior over time. The average speed (top) 
and fraction of time spent reversing (bottom) calculated over 100 s windows and 
averaged across individuals is shown for each strain. Shaded areas indicate the 
95% confidence interval. 
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Figure S5: Characterization of speed statistics for all strains. (A) The speed 
autocorrelation (black, calculated ignoring data within 1 s of reversal transi¬ 
tions) of each strain decays exponentially (red). (B) The speed distribution 
(black, calculated ignoring data within 1 s of reversal transitions) of each strain 
is closely reproduced by simulations of the model (red). Shaded areas indicate 
the 95% confidence interval. 
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Figure S6: Speed dynamics at reversal boundaries. (A) Cross-correlation 
function between the speed and reversal initiation or termination events (Sr) for 
individuals from the three example strains. Both the initiation and termination 
of reversals are associated with transient decreases in speed. In Figure]^ data 
within 1 s of these transitions were ignored (shaded gray region) and were not 
included as part of the model. (B) Comparison of a time constant for the 
reversal transition-speed cross-correlation, estimated from an exponential fit, 
with the Ts parameter for each trajectory (colored by strain). The black line 
indicates perfect correlation as expected if speed dynamics arose purely from 
reversal state transitions. Shaded areas indicate the 95% confidence interval. 


A 


$ JU775 QX1211 JU757 PS2045 DF5020 PS1159 sjh2 

¥ 5 



T (S) 


Figure S7: Orientation dynamics for all strains. (A) The mean absolute change 
in body orientation (black) increases linearly with time (red) in all strains. (B) 
The orientation autocorrelation (black) decays exponentially (red) in all strains. 
Shaded areas indicate the 95% confidence interval. 
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Figure S8: Analysis of orientation drift. (A-B) Analysis of orientation dynamics 
in an individual N2 worm. (A) The orientation over time is plotted along 
with the average drift {{k^{t))t). (B) The drift parameter k^{t) estimated in 
100 s windows is plotted over time. The measured drift changes sign over time. 
(C, top) The distribution of drift measurements (fc^) from B. (C, middle) The 
distribution of these measurements, averaged over time, across all N2 individuals 
is shown. The average drift is smaller and centered around 0. (C, bottom) To 
account for the larger orientation bias over 100 s windows, the absolute value of 
the bias in each window was taken and averaged over time. The distribution of 
this time-averaged magnitude of the bias is shown for all N2 individuals. 
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Figure S9: Analysis of reversals across all strains. (A) Distribution of A'lp 
for each strain shows two prominent peaks at 0 and 180 degrees. (B) The 
autocorrelation function of Atjj for each strain (black) along with the exponential 
fit (red). Shaded areas indicate the 95% confidence interval. 
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Figure SIO: Comparison of the average speed during forward and reverse runs. 
The average speed during forward and reverse runs is plotted for individual 
trajectories (dots, colored by strain). There is a strong correlation in the two 
speeds. 
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Figure Sll: Comparison of dimensionality reduction approaches. Two- 

dimensional projections of the behavior phenotypes of individual worms (circles) 
are shown for (A) PCA (as in the main text) and (B) t-stochastic neighbor em¬ 
bedding (t-SNE), a non-linear dimensionality reduction technique [^. Colors 
represent different strains. In (B), the long axis of the data is still well described 
by the top PCA mode. The gray dots are the predicted projections of the data 
points based on a linear fit of the top PCA mode to the 2D t-SNE projections. 
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Figure S12: Strain-specific phenotypic distributions. (A) The distribution of 
the average phenotype for each species {interspecies variation, blue) or individu¬ 
als within a strain (red) along the top behavioral mode is shown. The observed 
phenotypes used to estimate the distribution are indicated with dashes. (B) 
The bootstrapped standard deviation and 95% confidence interval of the dis¬ 
tributions in (A) are plotted. The average uncertainty in the species means or 
individual phenotypes is shown in the black bars. 
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Figure S13: Characteristics of behavioral phenotype distributions across 

strains. The distribution of projections of individual trajectories on the top 
behavioral mode for each strain was characterized by its mean and standard 
deviation. As a reference, the standard deviation of the average phenotype 
across species was also calculated. The mean and relative variability (standard 
deviation across individual trajectories normalized by the standard deviation 
across species) are plotted. C. elegans strains are shown in red. The error bars 
represent 95% confidence intervals. 
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Figure S14: Quantitative analysis of phenotypic divergence among strains. A 
metric distance for phenotypic distributions was calculated from the Jensen- 
Shannon distance 37 . The matrix shows the pairwise distance between the 


distributions of the mode 1 projections, shown in Figure Several clusters of 
similarity can be seen among related strains. For comparison, the phylogenetic 
tree is also shown. 



Figure S15: Overview of image processing steps. The video frames were pro¬ 
cessed by (1) subtracting the average of 50 frames evenly sampled from the entire 
movie and (2) cropping to each of the SDS-enclosed regions. (3) The largest 
worm-sized object was identihed following several image morphology operations, 
and (4) the centroid and image skeleton were measured. 
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Figure S16: Comparison of velocity calculation methods. Velocity autocorrela¬ 
tion functions for the three example strains with and without filtering of the data 
and without averaging over 100 s windows. The unfiltered velocity (black), esti¬ 
mated using a symmetrized derivative, contained a J-correlated short-timescale 
component in all strains that was particularly prominent in slow-moving strains 
such as PS312. The velocity calculated using a 1 s cubic polynomial hlter (red) 
does not contain this (^-correlated component. 
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Strain 

Mean 

(s) 

(fim/s) 

2.5% 

N2 

74 

62 

CB4856 

106 

88 

JU775 

no 

96 

QX1211 

46 

31 

JU757 

75 

57 

PS2045 

34 

25 

DF5020 

54 

36 

PS312 

24 

19 

PS1159 

30 

23 

sjh2 

162 

142 


.5% 

Mean 

DeS 

(fJLW? / s) 
2.5% 

97.5% 

88 

9800 

6100 

17100 

124 

26000 

13300 

39400 

124 

34700 

27200 

43000 

73 

1900 

500 

6300 

93 

2000 

800 

5500 

42 

5300 

1600 

11000 

69 

3400 

500 

9200 

28 

400 

200 

600 

41 

2000 

700 

5600 

186 

44100 

30200 

69600 


Table SI: Motility parameters. Each trajectory was fit to an effective random 
walk model described by speed ((s)) and an effective diffusivity (Z?eff)- Note that 
the speed parameter used here is lower than the model parameter due to the 
inclusion of reversal transitions. The distributions of these parameters were 
then bootstrapped to estimate the mean and the 2.5% and 97.5% percentiles 
(spanning the 95% confidence interval). The geometric rather than arithmetic 
mean was used as the parameters varied log-normally. 


Strain 

Mean 

ifim/s) 

2.5% 

97.5% 

Mean 

Ts 

(s) 

2.5% 

97.5% 

((M' 

Mean 

Ds 

mjs)'^ 

2.5% 

Is) 

97.5% 

N2 

94 

82 

105 

1.8 

1.4 

2.3 

500 

390 

650 

CB4856 

113 

96 

130 

1.5 

1.2 

2.0 

400 

320 

510 

JU775 

118 

102 

132 

1.4 

0.9 

1.9 

550 

360 

950 

QX1211 

74 

54 

107 

0.7 

0.6 

1.0 

810 

540 

1180 

JU757 

132 

107 

146 

2.2 

1.6 

2.9 

630 

440 

970 

PS2045 

64 

60 

69 

0.3 

0.3 

0.3 

1660 

1450 

1900 

DF5020 

78 

52 

93 

0.8 

0.6 

0.9 

1160 

490 

1800 

PS312 

36 

29 

42 

0.6 

0.4 

0.7 

550 

410 

740 

PS1159 

42 

31 

54 

0.6 

0.4 

1.1 

440 

240 

690 

sjh2 

170 

150 

194 

3.6 

2.9 

5.1 

760 

490 

1090 


Table S2: Speed dynamics parameters. The model parameters related to the 
speed dynamics are listed for each strain. For each worm in a strain, time- 
averaged parameters were calculated. This distribution was then bootstrapped 
to estimate the mean and the 2.5% and 97.5% percentiles of the distribution 
(spanning the 95% confidence interval). The geometric rather than arithmetic 
mean was used as the parameters varied log-normally. 
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Strain 

Mean 

1 

{rad/s) 
2.5% 

97.5% 

Mean 

Dip 

{rad? / s) 
2.5% 

97.5% 

N2 

0.039 

0.036 

0.042 

0.054 

0.047 

0.064 

CB4856 

0.044 

0.038 

0.051 

0.043 

0.034 

0.053 

JU775 

0.048 

0.041 

0.056 

0.045 

0.039 

0.052 

QX1211 

0.044 

0.036 

0.058 

0.055 

0.042 

0.079 

JU757 

0.039 

0.033 

0.045 

0.049 

0.040 

0.057 

PS2045 

0.038 

0.028 

0.058 

0.031 

0.021 

0.049 

DF5020 

0.034 

0.024 

0.046 

0.045 

0.033 

0.059 

PS312 

0.026 

0.022 

0.032 

0.044 

0.037 

0.056 

PS1159 

0.029 

0.021 

0.035 

0.026 

0.019 

0.033 

sjh2 

0.055 

0.049 

0.064 

0.118 

0.094 

0.143 


Table S3: Orientation dynamics parameters. The model parameters related 
to the orientation dynamics are listed for each strain. For each worm in a 
strain, time-averaged parameters were calculated. This distribution was then 
bootstrapped to estimate the mean and the 2.5% and 97.5% percentiles of the 
distribution (spanning the 95% confidence interval). The geometric rather than 
arithmetic mean was used as the parameters varied log-normally. 


Strain 

Mean 

"^fwid 

( 5 ) 

2.5% 

97.5% 

Mean 

"^rev 

( 5 ) 

2 . 5 % 

97 . 5 % 

N2 

13.2 

00 

bo 

20.3 

3.6 

2.9 

5.2 

CB4856 

29.6 

24.0 

39.4 

3.8 

2.8 

5.1 

JU775 

27.4 

22.3 

35.5 

3.3 

2.6 

4.2 

QX1211 

12.4 

4.9 

26.2 

1.9 

1.2 

3.4 

JU757 

4.5 

2.8 

8.2 

0.9 

0.7 

1.5 

PS2045 

7.3 

4.6 

18.1 

0.6 

0.4 

1.3 

DF5020 

6.3 

2.9 

10.3 

1.7 

0.9 

2.8 

PS312 

2.1 

1.4 

2.9 

1.5 

1.1 

2.0 

PS1159 

5.3 

3.2 

9.8 

0.9 

0.7 

1.2 

sjh2 

23.8 

18.5 

28.6 

6.5 

4.5 

9.6 


Table S4: Reversal state dynamics parameters. The model parameters related 
to the reversal state dynamics are listed for each strain. For each worm in a 
strain, time-averaged parameters were calculated. This distribution was then 
bootstrapped to estimate the mean and the 2.5% and 97.5% percentiles of the 
distribution (spanning the 95% confidence interval). The geometric rather than 
arithmetic mean was used as the parameters varied log-normally. 
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Parameter 

Mean 

Loading 

2.5% 

97.5% 

logio Ms 

0.27 

0.22 

0.32 

logio Ts 

0.28 

0.17 

0.39 

logio Ds 

-0.05 

-0.15 

0.04 

logio 141 

0.15 

0.12 

0.18 

logio 

0.17 

0.12 

0.24 

logio "^fwd 

0.72 

0.64 

0.77 

logio Trev 

0.52 

0.47 

0.58 


Table S5: Behavioral mode parameter loadings. The loadings of each parameter 
on the vector defining the top behavioral mode are listed. The values were 
bootstrapped as described in the main text to estimate the mean and the 2.5% 
and 97.5% percentiles of the distribution (spanning the 95% confidence interval). 


Strain 

Mean 

Projection 

2.5% 97.5% 

N2 

0.30 

0.06 

0.49 

CB4856 

0.56 

0.44 

0.70 

JU775 

0.50 

0.37 

0.62 

QX1211 

-0.01 

-0.44 

0.40 

JU757 

-0.30 

-0.55 

-0.02 

PS2045 

-0.62 

-0.92 

-0.21 

DF5020 

-0.28 

-0.72 

0.00 

PS312 

-0.78 

-1.01 

-0.59 

PS1159 

-0.60 

-0.93 

-0.23 

sjh2 

0.84 

0.69 

0.98 


Table S6: Behavioral mode projection by strain. The phenotypic projection 
along the first behavioral mode is listed for each strain. For each worm in a 
strain, a time-averaged projection was calculated. This distribution was then 
bootstrapped to estimate the mean and the 2.5% and 97.5% percentiles of the 
distribution (spanning the 95% confidence interval). 
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